#include "BlockedRadix.h"



void BlockedRadix::executeFFT() { executeFFT(false); }

void BlockedRadix::executeFFT(bool onlyComm) {


	int lgn = log2(numPoints);
	int lgp = log2(numNodes);

	int globalI = rank << (lgn - lgp);



	for (int s=0; s < lgn; s++) {

	
		// when s=0,numPoints=16 bitmask = 1    (1)
		// when s=3,numPoints=16 bitmask = 1111 (15)
		int bitmask =  (numPoints-1) >> (lgn-s-1);




		//distributed values
		if (s < lgp) {


			int partRank = rank ^ (1 << (lgp-s-1));
			MPI_Sendrecv(localData,2*numLocalPoints,MPI_DOUBLE,partRank,s,
				     localBuff,2*numLocalPoints,MPI_DOUBLE,partRank,s,
				     MPI_COMM_WORLD,&(comm->status));



		    	if (!onlyComm) {
				// U = localData, T = localBuff
				if (rank < partRank) {
			
					for (int i=0; i < numLocalPoints; i++) {
						int r = bitRev(i+globalI,numPoints);
						int omega_power = (bitmask&r) << ((lgn-1)-(s+1)+1);
						Complex omega = FFTEngine::omega(numPoints,omega_power);

						//update
						localData[i] += omega * localBuff[i];
					}
				} 

				// U = localBuff, T = localData
				else {
	
					for (int i=0; i < numLocalPoints; i++) {
						int r = bitRev(i+globalI,numPoints);
						int omega_power = (bitmask&r) << ((lgn-1)-(s+1)+1);
						Complex omega = FFTEngine::omega(numPoints,omega_power);

						//update
						localData[i] *= omega;
						localData[i] += localBuff[i];
					}
				
				}
			}
		}

		//local values
		else if (!onlyComm) {
			for (int i=0; i < numLocalPoints; i++) {
				int partner = i ^ (1 << (lgn-s-1));

				if (i < partner) {

					int r = bitRev(i+globalI,numPoints);
					//fprintf(stderr,"i=%d partner=%d r=%d r=%d\n",i,partner,r,r&bitmask);
					int omega_power = (bitmask&r) << ((lgn-1)-(s+1)+1);	
					Complex omega_t = FFTEngine::omega(numPoints,omega_power) * localData[partner];
					//fprintf(stderr,"omega_power=%d\n",omega_power);



					//update
					localData[partner] = localData[i] - omega_t;
					localData[i] = localData[i] + omega_t;	
				}
			}

		}
	}

	return;	
}

